In [1]:
import visual_behavior_glm
In [2]:
import visual_behavior_glm.src.GLM_params as glm_params
import visual_behavior_glm.src.GLM_analysis_tools as gat
from visual_behavior_glm.src.glm import GLM
import matplotlib.pyplot as plt
import visual_behavior.data_access.loading as loading
import visual_behavior.database as db
import plotly.express as px
import pandas as pd
import numpy as np
import os
from sklearn.cluster import KMeans
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
import seaborn as sns
In [3]:
%matplotlib inline
In [ ]:
get results for glm version 4¶
In [4]:
results_all = gat.retrieve_results(results_type='full')
In [5]:
results_all['glm_version'].unique()
#use v4
Out[5]:
In [6]:
rs = gat.retrieve_results(search_dict = {'glm_version': '4_L2_optimize_by_cell'}, results_type='summary')
len(rs)
Out[6]:
make a new column called 'identifier' that is {OPHYS_EXPERIMENTID}{CELL_SPECIMEN_ID}. This will give each a unique ID¶
In [7]:
rs['identifier'] = rs['ophys_experiment_id'].astype(str) + '_' + rs['cell_specimen_id'].astype(str)
In [8]:
rs
Out[8]:
In [ ]:
pivot the table to get it into a form for clustering¶
first get variance explained for full model to filter out poorly fit cells¶
In [9]:
model_output_type = 'variance_explained'
ve = rs.pivot(index='identifier',columns='dropout',values=model_output_type).reset_index()
ve
Out[9]:
In [10]:
cells_to_include = ve[ve['Full']>0.01].identifier.values
order = np.argsort(ve[ve.identifier.isin(cells_to_include)==True]['Full'])
cell_order = cells_to_include[order]
len(cells_to_include)
Out[10]:
get fraction change in explained variance after dropout of each regressor¶
In [11]:
model_output_type = 'fraction_change_from_full'
rsp = rs.pivot(index='identifier',columns='dropout',values=model_output_type).reset_index()
rsp
Out[11]:
In [12]:
tmp = ve.rename(columns={'Full':'varience_explained_full_model'})
rsp = rsp.merge(tmp[['identifier','varience_explained_full_model']], on=['identifier'])
filter out poorly fit cells¶
In [13]:
rsp = rsp[rsp.identifier.isin(cells_to_include)==True]
In [ ]:
merge in the identifier, cre_line, session_type and equipment_name columns¶
In [14]:
rspm = rsp.merge(rs[['identifier','cre_line','session_type','equipment_name']].drop_duplicates(),left_on='identifier',right_on='identifier',how='inner')
rspm
Out[14]:
add a column for session ID (1, 3, 4, or 6)¶
In [15]:
def map_session_types(session_type):
session_id = session_type[6:7]
return session_id
rspm['session_id'] = rspm['session_type'].map(lambda st:map_session_types(st))
rspm['session_id'].unique()
Out[15]:
save feature matrix to disk¶
In [16]:
# save = False
# if save:
# rspm.to_csv('/allen/programs/braintv/workgroups/nc-ophys/visual_behavior/ophys_glm/fraction_change_var_explained_v_4_L2_fixed_lambda=1_2020.08.09.csv', index=False)
get the column names to use for clustering¶
In [17]:
cols_for_clustering = [col for col in rspm.columns if col not in ['identifier','cre_line','session_type','equipment_name', 'session_id']]
cols_for_clustering = [col for col in cols_for_clustering if col not in ['image0','image1','image2','image3',
'image4','image5','image6','image7',
'visual']]
cols_for_clustering
Out[17]:
In [18]:
cols_for_clustering = [
'omissions',
'all-images',
'image_expectation',
'change',
'hits',
'misses',
'correct_rejects',
'false_alarms',
'post_lick_bouts',
'post_licks',
'pre_lick_bouts',
'pre_licks',
'rewards',
'pupil',
'running',
'time',
'model_bias',
'model_omissions1',
'model_task0',
'model_timing1D',
]
plot feature matrix¶
In [19]:
rspm[cols_for_clustering]
Out[19]:
In [20]:
feature_matrix = rspm[cols_for_clustering]
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(feature_matrix, vmin=-0.5, vmax=0.5, center=0, cmap='RdBu_r', ax=ax, cbar_kws={'label':model_output_type})
ax.set_ylabel('cells')
ax.set_title('GLM feature matrix')
Out[20]:
In [21]:
feature_matrix = rspm.sort_values('varience_explained_full_model').reset_index()[cols_for_clustering]
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(feature_matrix, vmin=-0.5, vmax=0.5, cmap='RdBu_r', center=0, ax=ax, cbar_kws={'label':model_output_type})
ax.set_ylabel('cells')
ax.set_title('GLM feature matrix')
Out[21]:
In [22]:
feature_matrix = rspm.sort_values('omissions').reset_index()[cols_for_clustering]
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(feature_matrix, vmin=-0.5, vmax=0.5, cmap='RdBu_r', ax=ax, cbar_kws={'label':model_output_type})
ax.set_ylabel('cells')
ax.set_title('GLM feature matrix')
Out[22]:
basic characterization¶
In [23]:
fig, ax = plt.subplots(figsize=(4,4))
ax = sns.scatterplot(data=rspm, x='omissions', y='all-images', hue='cre_line', ax=ax)
In [24]:
fig, ax = plt.subplots(figsize=(4,4))
ax = sns.scatterplot(data=rspm, x='omissions', y='pupil', hue='cre_line', ax=ax)
In [25]:
fig, ax = plt.subplots(figsize=(4,4))
ax = sns.scatterplot(data=rspm, x='pupil', y='running', hue='cre_line', ax=ax)
In [28]:
colors = sns.color_palette()
colors = [colors[0], colors[2], colors[3]]
cre_lines = np.sort(rspm.cre_line.unique())
In [35]:
for metric in cols_for_clustering:
fig, ax = plt.subplots(figsize=(6,4))
sns.pointplot(data=rspm, x='session_id', y=metric, hue='cre_line', hue_order=cre_lines, palette=colors, ax=ax)
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
PCA on the feature matrix¶
In [24]:
n_features = len(cols_for_clustering)
n_components = len(cols_for_clustering)
pca = PCA(n_components=n_components)
pca_result = pca.fit_transform(rspm[cols_for_clustering].values)
rspm['pc1'] = pca_result[:,0]
rspm['pc2'] = pca_result[:,1]
rspm['pc3'] = pca_result[:,2]
print('Explained variation per principal component: {}'.format(pca.explained_variance_ratio_))
plot explained variance for PCs¶
In [25]:
fig,ax=plt.subplots()
ax.plot(
np.arange(n_components),
pca.explained_variance_ratio_,
'o-k'
)
ax.set_xlabel('PC number')
ax.set_ylabel('variance explained')
ax.set_title('first 8 PCs explain >95% of the variance')
Out[25]:
In [26]:
np.cumsum(pca.explained_variance_ratio_)
Out[26]:
In [27]:
np.searchsorted(np.cumsum(pca.explained_variance_ratio_), .90)
Out[27]:
In [28]:
np.searchsorted(np.cumsum(pca.explained_variance_ratio_), .95)
Out[28]:
what are the PCs?¶
In [29]:
fig, ax = plt.subplots(figsize=(6,6))
ax = sns.heatmap(pca.components_, vmin=-1, vmax=1, cmap='RdBu_r', ax=ax, square=True,
robust=True, cbar_kws={"drawedges": False, "shrink": 0.7, "label": 'weight'})
ax.set_ylabel('principal components')
ax.set_xlabel('features')
# ax.set_title('principal axes in feature space \n(directions of maximum variance in the data)')
ax.set_ylim(0, n_components)
ax.set_xticklabels(cols_for_clustering, rotation=90);
In [30]:
pca.components_.shape
Out[30]:
In [31]:
fig,ax=plt.subplots(figsize=(12,4))
N_PCs = 8
for PC in range(N_PCs):
ax.plot(pca.components_[PC,:])
ax.legend(np.arange(N_PCs), title='PC', bbox_to_anchor=(1,1))
ax.axhspan(1/np.sqrt(32), -1/np.sqrt(32), zorder=-np.inf, alpha=0.25)
ax.set_xticks(np.arange(len(cols_for_clustering)))
ax.set_xticklabels(cols_for_clustering, rotation=45, ha='right')
ax.set_ylabel('weight')
ax.set_ylim(-1.1,1.1)
fig.tight_layout()
In [32]:
fig,ax=plt.subplots(figsize=(12,4))
for PC in range(8,20):
ax.plot(pca.components_[PC,:])
ax.legend(np.arange(10,21), title='PC', bbox_to_anchor=(1,1))
ax.axhspan(1/np.sqrt(32), -1/np.sqrt(32), zorder=-np.inf, alpha=0.25)
ax.set_xticks(np.arange(len(cols_for_clustering)))
ax.set_xticklabels(cols_for_clustering, rotation=45, ha='right')
ax.set_ylabel('??')
ax.set_ylim(-1.1,1.1)
fig.tight_layout()
plot covariance matrix - what does it mean? shouldnt omissions and images have negative covariance?¶
In [33]:
fig, ax = plt.subplots(figsize=(6,6))
ax = sns.heatmap(pca.get_covariance(), vmin=-0.2, vmax=0.2, cmap='RdBu_r', ax=ax, square=True,
robust=True, cbar_kws={"drawedges": False, "shrink": 0.7, "label": 'covariance'})
ax.set_title('covariance matrix')
ax.set_ylim(0, n_features)
ax.set_xticklabels(cols_for_clustering, rotation=90);
ax.set_yticklabels(cols_for_clustering, rotation=0);
plot PCA result - are these eigenvalues per cell? or original values * eigenvectors aka 'scores'?¶
In [34]:
pca_result[np.argsort(pca_result[:,0])]
Out[34]:
In [35]:
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(pca_result, cmap='RdBu_r', ax=ax,
robust=True, cbar_kws={"drawedges": False, "shrink": 0.7, "label": 'activation'})
ax.set_ylabel('cells')
ax.set_xlabel('PC')
ax.set_title('% variance explained in PC space')
ax.set_ylim(0, pca_result.shape[0])
ax.set_xlim(0, pca_result.shape[1])
ax.set_xticks(np.arange(0, pca_result.shape[1]));
# ax.set_xticklabels(cols_for_clustering, rotation=90);
In [36]:
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(pca_result[np.argsort(pca_result[:,0])], cmap='RdBu_r', ax=ax,
robust=True, cbar_kws={"drawedges": False, "shrink": 0.7, "label": 'activation'})
ax.set_ylabel('cells')
ax.set_xlabel('PC')
ax.set_title('% variance explained in PC space')
ax.set_ylim(0, pca_result.shape[0])
ax.set_xlim(0, pca_result.shape[1])
ax.set_xticks(np.arange(0, pca_result.shape[1]));
# ax.set_xticklabels(cols_for_clustering, rotation=90);
plot data in PC space¶
In [ ]:
In [37]:
# pc1 and pc2 columns in rspm correspond to pca_results for those PCs - is this the 'score' per cell?
fig,ax = plt.subplots(1,2,figsize=(12,6))
ax[0] = sns.scatterplot(data=rspm, x="pc1", y="pc2", hue="cre_line",
palette=sns.color_palette("hls", 3), legend="full", alpha=0.3, ax=ax[0])
# ax[0].set_xlim(-5,10)
# ax[0].set_ylim(-5,10)
ax[1] = sns.scatterplot(data=rspm, x="pc2", y="pc3", hue="cre_line",
palette=sns.color_palette("hls", 3), legend="full", alpha=0.3, ax=ax[1])
# ax[1].set_xlim(-5,10)
# ax[1].set_ylim(-5,10)
In [38]:
pca_result_df = pd.DataFrame(pca_result, index=rspm.identifier)
pca_result_df['cre_line'] = rspm['cre_line'].values
In [64]:
# pc1 and pc2 columns in rspm correspond to pca_results for those PCs - is this the 'score' per cell?
PC1 = 0
PC2 = 1
PC3 = 3
PC4 = 4
fig,ax = plt.subplots(1, 3, figsize=(15,5))
ax = ax.ravel()
i=0
ax[i] = sns.scatterplot(data=pca_result_df, x=PC1, y=PC2, hue="cre_line",
palette=sns.color_palette("hls", 3), legend="full", alpha=0.3, ax=ax[i])
# ax[i].set_xlim(-100,100)
# ax[i].set_ylim(-100,100)
i+=1
ax[i] = sns.scatterplot(data=pca_result_df, x=PC2, y=PC3, hue="cre_line",
palette=sns.color_palette("hls", 3), legend="full", alpha=0.3, ax=ax[i])
# ax[i].set_xlim(-100,100)
# ax[i].set_ylim(-100,100)
i+=1
ax[i] = sns.scatterplot(data=pca_result_df, x=PC3, y=PC4, hue="cre_line",
palette=sns.color_palette("hls", 3), legend="full", alpha=0.3, ax=ax[i])
# ax[i].set_xlim(-100,100)
# ax[i].set_ylim(-100,100)
fig.tight_layout()
In [ ]:
In [41]:
query_string = '''pc1>-100 and pc1<100 and pc2>-100 and pc2<100 and pc3>-100 and pc3<100'''
fig = px.scatter_3d(
rspm.query(query_string),
x='pc1',
y='pc2',
z='pc3',
color='cre_line',
)
fig.update_traces(
marker=dict(
size=3,
opacity=0.25
)
)
fig.update_layout(
margin=dict(l=30, r=30, t=10, b=10),
width=1200,
height=1000,
)
# fig.write_html("/home/dougo/code/dougollerenshaw.github.io/figures_to_share/2020.08.09_PCA_on_GLM.html")
# fig.show()
In [ ]:
In [42]:
# # pc1 and pc2 columns in rspm correspond to pca_results for those PCs - is this the 'score' per cell?
# fig,ax = plt.subplots(n_components, n_components, figsize=(20,20))
# ax = ax.ravel()
# i = 0
# for PC1 in range(n_components):
# for PC2 in range(n_components):
# ax[i] = sns.scatterplot(data=pca_result_df, x=PC1, y=PC2, hue="cre_line",
# palette=sns.color_palette("hls", 3), legend="full", alpha=0.3, ax=ax[i])
# ax[i].set_xlim(-100,100)
# ax[i].set_ylim(-100,100)
# # ax[1] = sns.scatterplot(data=rspm, x="pc2", y="pc3", hue="cre_line",
# # palette=sns.color_palette("hls", 3), legend="full", alpha=0.3, ax=ax[1])
# # ax[1].set_xlim(-100,100)
# # ax[1].set_ylim(-100,100)
In [ ]:
cluster with k-means¶
cluster on feature matrix¶
In [43]:
feature_matrix = rspm[cols_for_clustering]
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(feature_matrix, vmin=-0.5, vmax=0.5, cmap='RdBu_r', ax=ax, cbar_kws={'label':model_output_type})
ax.set_ylabel('cells')
ax.set_title('GLM feature matrix')
Out[43]:
In [59]:
fig, ax = plt.subplots(figsize=(4,4))
ax = sns.scatterplot(data=rspm, x='omissions', y='all-images', hue='cre_line', ax=ax)
# ax.set_xlim(-100, 100)
# ax.set_ylim(-100, 100)
In [58]:
fig, ax = plt.subplots(figsize=(4,4))
ax = sns.scatterplot(data=rspm, x='omissions', y='running', hue='cre_line', ax=ax)
# ax.set_xlim(-100, 100)
# ax.set_ylim(-100, 100)
In [57]:
fig, ax = plt.subplots(figsize=(4,4))
ax = sns.scatterplot(data=rspm, x='omissions', y='pupil', hue='cre_line', ax=ax)
# ax.set_xlim(-100, 100)
# ax.set_ylim(-100, 100)
In [47]:
kmeans = KMeans(n_clusters=10)
kmeans_result = kmeans.fit_predict(feature_matrix)
rspm['kmeans_result_features'] = kmeans_result
rspm['kmeans_result_features'].value_counts()
Out[47]:
In [ ]:
cluster on PCA results¶
In [48]:
fig, ax = plt.subplots(figsize=(6,10))
ax = sns.heatmap(pca_result, cmap='RdBu_r', ax=ax,
robust=True, cbar_kws={"drawedges": False, "shrink": 0.7, "label": '?'})
ax.set_ylabel('cells')
ax.set_xlabel('PC')
ax.set_title('')
ax.set_ylim(0, pca_result.shape[0])
ax.set_xlim(0, pca_result.shape[1])
ax.set_xticks(np.arange(0, pca_result.shape[1]));
# ax.set_xticklabels(cols_for_clustering, rotation=90);
In [49]:
kmeans = KMeans(n_clusters=10)
kmeans_result = kmeans.fit_predict(pca_result)
rspm['kmeans_result'] = kmeans_result
rspm['kmeans_result'].value_counts()
Out[49]:
"Then we applied consensus clustering to the PCs, by running K-means using the PCs 100 times until reaching a stable co-clustering association matrix, where each entry represents the probability of two units belonging to the same cluster." - Xiaoxuan paper
clustering on the PCs - is this legit?¶
In [50]:
fig, ax = plt.subplots(figsize=(6,6))
ax = sns.heatmap(pca.components_, vmin=-1, vmax=1, cmap='RdBu_r', ax=ax, square=True,
robust=True, cbar_kws={"drawedges": False, "shrink": 0.7, "label": 'contribution'})
ax.set_ylabel('principal components')
ax.set_xlabel('GLM features')
ax.set_title('principal axes in feature space \n(directions of maximum variance in the data)')
ax.set_ylim(0, n_components)
ax.set_xticklabels(cols_for_clustering, rotation=90);
In [51]:
kmeans_result
Out[51]:
In [52]:
kmeans = KMeans(n_clusters=10)
kmeans_result = kmeans.fit_predict(pca.components_)
kmeans_result
Out[52]:
In [ ]:
Agglomerative clustering¶
In [53]:
# import hierarchical clustering libraries
import scipy.cluster.hierarchy as sch
from sklearn.cluster import AgglomerativeClustering
on features¶
In [54]:
# create clusters
hc = AgglomerativeClustering(n_clusters=4, affinity = 'euclidean', linkage = 'ward')
# save clusters for chart
y_hc = hc.fit_predict(rspm[cols_for_clustering])
In [55]:
rspm['hc'] = y_hc
rspm['hc'].value_counts()
Out[55]:
on PCA results¶
In [56]:
# create clusters
hc = AgglomerativeClustering(n_clusters=4, affinity = 'euclidean', linkage = 'ward')
# save clusters for chart
y_hc = hc.fit_predict(pca_results)
In [ ]:
rspm['hc_pca'] = y_hc
rspm['hc_pca'].value_counts()
In [ ]:
In [ ]: